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Abstract 

In  this  paper  we  study  the  transient  behavior  of  the  A/GEi/A/GE.vf /I 
queueing  system,  where  MGE  is  the  class  of  mixed  generalized  Erlang  distri- 
butions which  can  approximate  an  arbitrary  distribution.  We  use  the  method 
of  stages  combined  with  the  separation  of  variables  and  root  tinding  techniques 
together  with  linear  and  tensor  algebra.  We  find  simple  -iDsed  form  expres- 
sions for  the  Laplace  transforms  of  the  queue  length  distribution  and  the  wait- 
ing time  distribution  under  FCFS  when  the  system  is  initially  empty  and  the 
busy  period  distribution.  \Ve  report  computational  results  by  inverting  these 
expressions  numerically  in  the  time  domain.  Because  of  the  simplicity  of  the 
expressions  derived  our  algorithm  is  very  fast  and  robust. 
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1      Introduction 

Transient  and  busy  period  analysis  in  queueing  models  have  long  been  considered 
as  very  difficult  problems.  Yet,  in  many  situations  it  is  very  important  to  study 
the  transient  behavior  of  queueing  systems.  For  example,  systems  often  encounter 
transient  behavior  due  to  exogenous  changes,  such  as  the  opening  or  closing  of  a 
queueing  system  or  the  application  of  a  new  control.  Furthermore,  even  in  systems 
with  time  homogeneous  behavior  the  convergence  to  steady  state  is  so  slow  that  the 
equilibrium  behavior  is  not  indicative  of  system  behavior.  Examples  from  practical 
situations,  in  which  transient  phenomena  are  important,  include  manufacturing  sys- 
tems with  frequent  start  up  periods  and  transportation  systems  with  time  varying 
demand  (for  example  airport  runway  operations  in  major  airports). 

Analytical  investigations  of  the  transient  behavior  of  queueing  systems  are  very 
rare,  mainly  because  of  the  complexity  involved.  For  the  M/M/\  queue  expressions 
for  the  queue  length  probabilities  are  known  as  sums  of  modified  Bessel  functions 
(see  Gross  and  Harris  [6]).  An  indication  of  the  interest  of  the  research  community 
in  transient  behavior  of  queues  is  the  recent  work  of  Abate  and  Whitt  [1]  for  the 
M/M/1  queue. 

In  the  last  decade,  work  on  the  transient  behavior  of  queueing  systems  has 
concentrated  on  numerical  techniques.  This  change  in  emphasis  was  primarily  mo- 
tivated by  the  analytical  complexity  of  the  problems  involved.  The  two  principal 
methods  are  the  randomization  technique  introduced  by  Grassmann  [5]  and  numer- 
ical integration  methods  of  the  underlying  Kolmogovov  differential  equations  (see 
Gross  and  Harris  [6],  section  7.3.2  and  references  therein). 

In  this  paper  we  study  various  transient  performance  characteristics  of  the 
MGEi,/MCE\f/l  system,  which  is  an  important  special  case  of  the  CI/G/l  queue. 
In  a  sequel  paper  (Bertsimas  et.  al.  [3])  we  formulate  the  problem  of  finding  si- 
multaneously the  waiting  time  distribution  and  the  busy  period  distribution  of  the 
GI/G/l  queue  with  arbitrary  distributions  as  a  Hilbert  problem.  The  MGE,  which 


is  described  in  some  detail  in  Section  2,  is  the  class  of  mixed  generalized  Erlang 
distributions,  which  is  dense  in  the  space  of  all  distributions  and  can  approximate 
arbitrarily  closely  every  distribution.  For  a  discussion  of  the  properties  of  the  MGE 
class  see  Bertsimas  [2].  We  study  the  queue  length,  the  waiting  time  and  the  busy 
period  distribution.  We  use  the  separation  of  variables  technique  together  with  root 
finding  techniques  to  establish  closed  form  expressions  for  the  Laplace  transform 
of  the  distributions  under  study  (queue  length,  waiting  time,  busy  period).  The 
advantage  of  these  closed  form  expressions  is  that  they  are  relatively  simple  and  can 
be  used  for  numerically  inverting  them  in  the  time  domain.  In  fact,  in  Section  6 
we  report  computational  results  for  inverting  numerically  the  transform  of  the  busy 
period  distribution. 

These  expressions  also  explain  the  difficulty  that  the  research  community  has 
had  over  the  years  in  establishing  expressions  for  the  distributions  we  study  in  the 
time  domain.  Despite  their  simplicity  in  the  transform  domain,  our  expressions 
involve  roots  of  polynomial  equations.  In  general,  these  roots  can  not  be  computed 
analytically  and  even  if  they  are  known  they  are  complicated  enough  to  make  their 
analytic  inversion  extremely  complicated  if  not  impossible. 

The  paper  is  structured  as  follows.  In  Section  2  we  describe  the  MGE  distri- 
bution and  the  notation  we  use.  In  Section  3  we  derive  closed  form  expressions  for 
the  transform  of  the  queue  length  distributions  when  the  system  is  initially  empty. 
In  Section  4  we  find  an  explicit  expression  for  the  busy  period  distribution  while 
in  Section  5  we  analyze  the  waiting  time  distribution  under  FCFS.  In  Section  6  we 
describe  the  algorithm  to  invert  numerically  the  closed  form  expressions  derived  in 
the  previous  sections  and  also  report  some  preliminary  computational  results.  The 
final  section  contains  some  concluding  remarks. 


2      Notation 

The  general  Coxian  class  Cn  was  introduced  in  Cox's  [4]  pioneering  paper.  The 
stage  representation  of  the  Coxian  distribution  is  presented  in  Figure  1.  It  should 
be  noted  that  this  stage  representation  of  the  Coxian  distribution  is  purely  formal 
in  the  sense  that  the  branching  probabilities  q,  can  be  negative  and  the  rates  /x, 
can  be  complex  numbers.  The  mixed  generalized  Erlang  distribution  is  a  Coxian 
distribution,  where  we  assume  that  the  probabilities  7,  are  non-negative  and  the 
rates  n,  are  reals.  As  a  result,  the  mixed  generalized  Erlang  distribution  has  a  valid 
probabilistic  interpretation,  which  is  further  exploited  in  this  paper. 
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Figure  1:  The  MCE\f  class  of  distributions 

To  analyse  the  model  we  conceive  of  the  arrival  process  as  an  arrival  timing 
channel  (ATC)  consisting  of  L  consecutive  stages  with  rates  Aj ,  A2, ....  A/,  and  with 
probabilities  Pi,P2,  -^PL  =  1  of  entering  the  system  after  the  completion  of  the 
1st,  '2nd,  ,  .  ,  Lth  stage.  We  remark  that  as  soon  as  a  customer  in  the  ATC  enters 
the  system  a  new  customer  arrives  at  stage  1  of  the  ATC,  For  the  service  time 
distribution  we  consider  as  above  a  service-timing  channel  (STC)  consisting  of  M 
consecutive  stages  with  rates  Mi^ /^2i  •1  A*M  ^"d  with  probabilities  <?i,  92.  ■••.  9A/  =  1 
of  leaving  the  system. 

In  the  transient  domain  we  introduce  the  random  variables: 


N{t)  =  The  number  of  customers  in  the  system  at  time  t. 


Ra{t)  =  The  ATC  stage  currently  occupied  by  the  arriving  customer  at  time  t. 

R,[t)  =  The  STC  stage  currently  occupied  by  the  customer  who  is  being  served  at 
time  t. 

W{t)  —  The  waiting  time  of  a  customer  arriving  at  time  t. 

With  the  above  definitions  the  system  can  be  formulated  as  a  continuous  time 
Markov  chain  with  infinite  state  space: 

{(A'(f),/2a(t).«5(O)..V(t)  =  0.  1 fla(n=  1.2 L.    /?,(<)=  1.2 A/}. 

We  now  introduce  the  following  set  of  probabilities: 

Pn,..,(f)  =  Pr{.V(i)  =  n.R,{\)  =  i.R,{i)^j\.       P^Jt)  =  Pr{.\(t)  =  0.  Ra{t)  =  i} . 
We  will  also  use  the  following  notation: 

•  T  =  the  mean  interarrival  time.  -  =  the  mean  service  time  and  p  =  -  =  the 
traffic  intensity. 

•  P^,{t)  =  is  a  column  vector,  whose  elements  are  the  probabilities  Pn.i.j[t)- 

•  ak{t).br{i)=the  probability  density  function  (pdf)  of  the  remaining  interarrival 

(service)  time  if  the  customer  in  the  ATC  (STC)  is  in  stage  k   =   1 L 

(r  =  1,.....\/).    Therefore,  ai(t).6i(t)  is  the  pdf  of  the  interarrival  (service) 
time. 

•  a^(i).tf(t)=  the  probability  to  move  from  stage  i  <  j  of  the  ATC  (STC)  to 
stage  j  during  the  interval  t  without  having  any  new  arrival  (service  comple- 
tion). 

.  a-,{t)  =  {a\{t),...,a^{t)Y,a-l{t)  =  (0 at(t), . . .  ,a^(t))'. 

•  h{t)  =  {b\{t),....blf{t)y.b1{t)  =  {0,....blit),...,bl'{t)y, 

•  Tfn{s),ak{s),3k{s)  =  the  Laplace  transforms  of  P„(i),  aX(t)i  ^fc(0  respectively. 


.  e- =(o,...,i,...,oy 

By  introducing  the  following  upper  semidiagonal  matrices  Aq,    Bq  and  the  dyadic 

matrices  Ai ,   fii: 

*"  A,     -(l-p,)>i     •••      0 
.4o  = 


.4,  = 


-PiAi      0 


-Az,       0      ■•■    0 
and  similarly  Bq,  B\  for  the  service  time  distribution  we  have 


Q-;'(s)  =  e-;(75  +  .4o)-', 

L  L 


n::ai-p.)A. 


Qfc(s)  =  -ei.{Is+  An)    '.4,^1  =  S^prKalis)  -  Y^PrAr     '/    ,         ' 

and  similarly  for  (3^  (s),  dk{s)- 

Finally  we  use  the  usual  tensor  notation  (see  also  Neuts  [10],  p. 53) 

(ll,...,In)'®  (yi,..  .  ,ym)'  =  (-Ciyi,--  •■•Tnyi-  ■■■,Xr^ym)'  ■ 

Finally,  ii"'",£°,£~  .^,£7,^  are  transition  rate  matrices.  We  can  express  the  tran- 
sition rate  matrices  in  terms  of  the  matrices  Aq,  A\,  Bq,  B\  using  tensor  notation. 
For  example 

£°  =  -Aq®  I  -  I®  Bq. 

3      The  queue  length  distribution  in  the  transient  do- 
main 

The  system  MCEi/MGE\i/\  is  an  instance  of  the  homogeneous  row-continuous 
Markov  chain   with  a  single  boundary  (Keilson  and  Zachmann  [8]).     We  analyze 
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the  homogeneous  part  of  the  Markovian  dynamics  using  the  separation  of  variables 
technique  combined  with  tensor  algebra.  Then  we  analyze  the  compensation  part 
(i.e.  the  boundary  condition  plus  the  initial  condition)  using  linear  algebra.  In 
this  way  we  succeed  in  finding  a  closed  form  expression  of  the  Laplace  transform 
of  the  queue  length  distribution.  We  assume  that  the  system  is  initially  empty. 
Although  our  approach  can  be  in  principle  applied  even  in  the  case  of  arbitrary 
initial  conditions,  the  algebra  required  makes  the  explicit  derivation  very  hard. 

The  homogeneous  part 

Using  the  notation  of  section  2  we  first  write  the  Chapman-Kolmogorov  forward 
equation  for  n  >  1: 

By  taking  the  Laplace  transform  and  using  the  assumption  that  Pn{0)  =  0  (for 
n  >  1)  we  obtain: 

s?;(s)  =  ir;_i(sk"^  +  ^n(s)g°  +  ^n+iis)g~  (1) 

The  compensation  part 

Similarly,  the  Chapman-Kolmogorov  forward  equation  for  n  =  0  and  n  =  1  is: 

J^Poi^)  =  Poit)^  +  Pirn: 

and  therefore  the  Laplace  transform  of  the  the  equations  for  n  =  0,  1  is: 

s^o{s)  =  Po{0)  +  ^'o{s)^  +  ^[{s)g;  (2) 

5?i(s)  =   ro(s)^  +  ^xis)g°  +  ^2{s)g-  (3) 

Our  general  strategy  for  analyzing  these  equations  is  the  following.    We  first 
find  the  general  solution  of  (1)  for  n  >  1.    The  solution  of  the  Laplace  transform 


Tfn(s)  is  then  written  as  a  linear  combination  of  M  geometric  terms.     The  only 

unknowns  are  the  M  constants  that  depend  on  s  (the  coefficients  Dr,r  =  [ M 

below).  From  (2)  we  find  Tfo{s)  as  a  function  of  the  M  coefficients  D^.  Finally 
we  use  (3)  to  find  a  linear  system  of  M  equations  with  A/  unknowns.  W'e  exploit 
the  particular  structure  of  the  linear  system  to  find  a  closed  form  solution  for  the 

unknowns  Dr,  r  =   1 M     .\s  a  result,  we  find  explicit  closed  form  expressions 

tor  the  Laplace  transform  of  the  queue  length  distribution. 

The  advantage  of  our  expressions  is  that  they  can  be  numerically  inverted  in 
real  time.  In  principle  this  approach  works  for  arbitrary  initial  conditions.  We  were 
able  to  find  closed  form  expressions  only  in  the  case  in  which  the  system  is  initially 
empty,  i.e.  P„(0)  =  0  for  n  >  1.  We  can  now  prove  our  first  result. 

Theorem  1  Under  the  assumption  that  the  system  is  initially  empty,  i.e.  the  only 
nonzero  initial  probabilities  are  Po,fc(0),  and  p  <  1  the  transform  of  the  queue  length 
distribution  has  the  following  form: 

M 

wn{s)     =     ^DrPi{xr{s))®ai{s-Xr{s))iai{s-Xr{s)))^-^   {n>l} 

L  M  , 

Ms)     =     ^Po.fc(0)a";(5)  +  ^D,— -,Ji(x,(s))(Q-i(5-j:,(s))-ai(s)), 
k  =  l  r=l         ^'-^^) 

where 

p    _i:LxPo.kiO)<yk{s){-l)''3^'{0)  it  x,{s) 

!-«,(.)  ;3'^'{xr{s))    ""^^'^y  xr{s)-x,{s)' 

and  X  =  Xr{s)  (for  r  =  1,  ..  .  ,  MJ  are  the  .\f  roofs  of  the  equation: 

r  ai(s-x)/?i(x)=  1 

l3fJ(x)  <  0(i.e.  ||ai(s- j-)||  <  1)   for3?(s)  >  0 

Proof 

We  first  find  the  general  solution  of  (I).  Equations  (1)  are  partial  difference  equa- 
tions with  constant  coefficients.  Following  the  separation  of  variables  technique  we 
assume  that  a  general  solution  Tn(s)  of  the  difference  equation  (1)  is: 

;r~(s)  =  v..(s)t^,(s)(ti;(s))"-\ 
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which  can  be  written  in  tensor  notation  as  follows: 


T„(S)  =   ^{S)  ®  11'{S){W{S)) 


n-1 


With 


^{s)  = 


\   rL{s)  I 


Ais)  = 


Vi(*)   \ 


\   Vm[s)  I 


We  substitute  the  form  of  ?n(5)  into  (1)  and  using  tensor  notation  we  obtain: 

s^(s)  ®  ^'{s){w{s)f-^     =  -/(5).4i  ®  v^'(s)(u(s))"-2  -  /(5)^o  ®  t/i''(s)(a.(s)r-' 

-  ^{s)  ®  ^'(s)Bo(it<s)r-i  -  /(s)  ®  ^~(s)5,(u.•(5))^ 

which  by  collecting  terms  can  be  written  as 

s(p{s)  ®  ^'{s)  +  ^{s){Ao  +  -)-.4i )  ®  t/r'(s)  +  ^(s)  ®  .?(s)(flo  +  u.'(s)fli)  =  0. 

Using  the  standard  separation  of  variables  arguments  (see  Bertsimas  [2])  in  tensor 
notation  we  demand  that  ^(s)  is  a  row  eigenvector  of  the  matrix  (.4o  +  ^^TTT  -^i) 
with  an  eigenvalue  —y{s)  and  v'l*)  is  a  row  eigenvector  of  the  matrix  [Bq  +  wi^s)B\) 
with  an  eigenvalue  — x^s).  As  a  result, 

(s-x(s)-y(s))[^(s)®T^'(s)]  =  0, 


and  therefore 


s  =  i{s)  +  y{s). 


(4) 


In  the  following  claim  we  establish  the  relation  among  Ui{s).  y(s)  and  x{s). 
by  computing  the  characteristic  polynomials  of  the  matrices  (.4o  -i-  ^7777 -^i)  and 
(Bo  +  u(s)fli)- 

Claim  1   Qi{t/(s))    =    w{s)  and  w{s)(3\{x{s))    =    1. 


Proof 


Since  x{s)  is  an  eigenvalue  of  (flo  +  w{s)Bi),  it  satisfies  the  following 
characteristic  equation: 

det[/r(s)  +  flo  +  ti'(s)Si]  =  0, 

that  is 

det[/r(s)  +  Bo]  det[/  +  w{s){Ix{s)  +  flo)"'5i]  =  0. 

Since  the  matrLV  Bq  has  full  rank,  det(/ir(s)  +  flo)  =  Denominator [Ji(x(s))]  = 
UrLiif^r+i{s))  ^  0.  Also  for  any  rank  1  matrix  S,  det(/+S)  =  1  +  trace(S). 
Thus,  since  Bi  has  rank  1 

1  +  u(s)trace((/r(s)  +  Bo)~'fii  =  0. 

But  /?i(i(s))  =  -  trace((/r(s)  +  Bq)-^ Bi)  and  therefore 

u(5)A(r(s))=l. 

Using  exactly  the  same  methodology  we  establish  that 

ai(y(s))  =  u-(s).a 

We  now  compute  in  the  following  claim  the  eigenvectors  li'(s)  and  ifi{s)  in  closed 
form. 

Claim  2  The  vector  cri(y(s))  is  a  row  eigenvector  of  the  matrix  [Aq  +  ^^TTf  -"^i) 
and  the  vector  3i{x(s))  is  a  row  eigenvector  of  the  matnx  (  Bq  +  u(s)Si).  Therefore 
we  can  choose  'fi{s}    =    0*1(1/(5))  and  xp{s)    =    /?i(x(s)). 

Proof 

We  prove  the  claim  for  i3i{x{s)).    The  case  for  0*1(1/(5))  is  similar. 
Since  i*i  {x{s))  =  e*,(/j(s)  +  5o)~'  we  have  that 

/fl'(i(s))(fio  +  u'(5)5i)  =     e-i(/j(5)  +  BorHBo  +  u(5)Si) 

=     e-'i(/ir(5)  +  Bo)-\-Ix{s)  +  (/x(s)  +  ^o)  +  a'(5)5i) 
=     -xis)i3iix{s))  +  e-'i  -  uis)3,{xis))e'„ 
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since  Jj  (s)fli   =  — i3i(s)e\.    But  from  claim  1,  u,'(s)5i(i(s))  =   1  and 
thus 

l3i{x{s)){Bo  +  w{s)Bi)  =  -x{s)3i\x{s)).n 

As  a  result  of  (4)  and  claim  1,  and  since  we  are  looking  for  roots  w{s)  =  ai(y(s))  = 
Qi(s  —  r{s))  inside  the  unit  circle,  so  that  the  solution  is  stationary,  the  general 
solution  of  equation  (1)  is 

^n{s)  =  3,{x{s))  S  a,{s  -  x{s))  {a,{s  -  x{s))f-'  ,  (5) 

where,  because  of  (4)  and  claim  I.  x  =  i[s)  satisfy  the  equations 


(6) 


Qi(s  -  x)3i{x)  =  1 
3?(r)  <  0  (i.e.  |iai(s  -  r)j|  <  1)   for  ^{s)  >  0 

In  the  following  claim  we  investigate  the  number  of  roots  of  (6). 

Claim  3    For  p  <  1  equation  (6)  has  M  roots  Xr{s).   r  =  1 A/. 

Proof 

This  can  be  easily  established  from  an  application  of  Rouche's  theo- 
rem in  the  domain  ^{x)  <  0.  For  a  very  similar  application  of  Rouche's 
theorem  see  Bertsimas  [2].  The  same  result  can  be  established  from  ma- 
trix geometric  considerations  by  noticing  that  the  roots  ai(s— Xr(s)),  r  = 
1, . . . ,  A/  are  the  A/  eigenvalues  of  the  matrix  R{s)  in  Ramaswami  [12]. □ 

Assuming  that  the  A/  roots  of  (6)  are  distinct  we  can  now  write  an  explicit  expression 
for  T„(s),n  >  1  by  taking  linear  combinations  of  the  general  solution  form.  Indeed, 
there  are  coefficients  Dr.r  —  1 \I  such  that 

\t 

f„(s)  =  Y.  Dji{Xr{s))  ®  ai{s  -  Xr{s))  (qi(s  -  Xris})^-'   {n  >  1} .  (7) 

r  =  l 

Remark:   The  distinctness  assumption 

This  distinctness  assumption,  is  very  convenient  in  order  to  find  an  explicit  expres- 
sion for  7rn(s),  n  >  1.  but  it  is  not  critical  however.  The  algebraic  theory  of  rational 
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functions  guarantees  that  if  there  are  multiple  roots,  we  can  take  the  limit  of  (7)  as 
Xr{s)  — '  ikis)  for  some  r.k.  In  other  words,  we  first  solve  the  problem  assuming 
that  the  roots  are  distinct  and  at  the  final  stage  we  show  the  results  are  independent 
from  this  assumption.  In  fact,  our  final  expressions  for  the  queue  length  and  the 
busy  period  distributions  are  simple  symmetric  functions  of  these  roots.  So,  finding 
the  limit  in  the  case  where  there  are  multiple  roots  is  an  eeisy  task. 

The  remaining  unknowns  are  the  coefficients  D^  (r  =  1,...,M)  and  tto{s).    In 
the  following  claim  we  express  Tro{s)  as  a  linear  combination  of  the  Dr's. 

Claim  4 

L  M  . 

?o(s)  =  Y.  PoA^Wkis)  +  YL  Dr^^M^ris))  (ai(s  -  x,(s))  -  ai(s))  .        (8) 

fc=I  r=l  ■^'•V*^ 

Proof 

We  substitute  (7)  into  the  equation  (2)  and  we  obtain 

A/ 

SK'ais)  =  P^(0)  -  ro(s).4o  -  J]  a  (Ji(x.(s))Si.-,)  ai'(5  -  x.(s)). 

r=l 

Thus 

M 

^0(5)     =     PoWils  +  Ao)-'  +  J2  Dr3,{xr{s))e',ilis  -  xris))  +  Ao}-\ls  +  Ao)-' 

L  M  . 

=       Y  /'O,)c(0)a-;(S)  +  J2  Dr—--f3i{Xris))  (ori(s  -  Xr{s))  -  0.1(5))  .D 
k=l  r  =  l  ^^-y^' 

The  next  step  in  our  approach  is  to  find  the  coeflRcients  Dr  (r  —  I,...,A/).  In 
the  following  claim  we  establish  the  equations  from  which  the  coefficients  Dr  are 
computed. 

Claim  5    For  all  k  =  I M: 

M 

Proof 
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y^^      1      l3^{xAs))  ^      T.k=i  Po.>c{0)c^k{s)  ^g^ 


Using  (7)  we  easily  obtain  that 

Af 

S7r[{s)  -  /i(s)£°  -  T2(s)^~    =  Y^  Drth   (jr(s))- 

r=:l 

From  (3) 

M 

Y^  DrMxris))  =  -(/o(s)^ie-j)e-i.  (10) 

r  =  l 

Using  (8)  and  since  q"!  (s)(  — .4ie*i)  =  Q\{s)  we  obtain  that 

'^'  1 

-(T'o(s).4ie-,)     =     Po(0)(/s  +  Ao)-'Aiei  +  T  D,-^/5i(x,(s))  {a,{s  -  x.(s))  -  a,{s)) 

where  we  defined  a'{s)  =  Po(0)(/s  +  ^o)~'-  We  substitute  this  equation 
in  (10)  and  obtain 

'r   \-<  STniT't      (^^        1  -  Qi  (s)^i  (Jr(5))  ^  , 

=  x:aj,(x.(.)){/---^-^-«,(.)-^}. 

Since  Ji'(s)  =  e*i(/s+  So)"'  and  Ji(s)e*i  =  -;i*i'(s)5i,  then 

.v/ 

^  n  

.(s) 

1      -. 


c''is)e''i     =     TD,—--J^'{xr{s)){Ixr{s)-{IXr{s)  +  Bo)-a^is)Bi} 

M 

=     -YDr-^/i'(xr{s)){Bo  +  ax{s)Bi). 


But  e"\(/s  +  5o  +  ir5i)-i  =  ,/t'^4'^(,)  and  hence  €"\(5o  +  ai(s)Si  )-^  = 
'  \' \  .  As  a  result, 

Yo'A'ix.is))  =  -^llL.^r(0) 

;fr;      •rr(s)  i-ai(s) 

=   -%iM5)^/,'(o),      (U) 

l-ai(s) 
Therefore,  for  all  k   —    1,...,A/  the  coefficients  Dr  satisfy  the  linear 
system  (9).  □ 
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We  now  solve  the  linear  system  (9)  in  closed  form. 


Claim  6    The  linear  system  (9)  has  the  following  solution:  For  all  r    —    1, . . . ,  A/.- 

^    i:LxP^,mc'k(s){-\)''3^'{^)  .  .^     x,{s) 


l-a,(.)  JF{xAs)) 


k=l 


Xr{s)  -  Xfc(s) 


Proof 


Let  Cr  =  ^L   '""'^'^ 4fe^D,.  Since  for  all  i    =    1,...,A/ 

we  obtam  that  for  all  t  =  1 A/,  J^^^.  Cr  ow'^'^n  =  1.  '•«■ 


r  =  l  u  =  fc  +  l 


^^n 


We  expand  the  above  equation  as  a  polynomial  of  rr{s)  and  obtain 


J2    E    Cr<T,A-^r{s)r   =   1, 


r  =  l   n=0 

where  (j^.n  are  the  coefficients  in  the  expansion.  We  express  this  equation 
in  matrix  form  and  we  obtain 


1 


"■i.o 


(^i(*^)) 


w-i 


{xm{s)) 


Af-1 


0     ■••       1 

s 

Since  by  the  definition  of  the  (7^. 


a- 

A 


(  c,  \ 

( l\ 

\  ^'^' ) 

[  1 1 

Jl            u 

c 

— 

1 

'^'l.O 


/  r\f-i  \        /  n^'-i 


n^LV(l  +  ;i;^)\ 


1      V    I    /     V 
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By  letting  x  =  0.  we  obtain  that  S     1  =  e\i  and  thus, 

We  now  observe  that  the  matrix  A  is  a  Vandermode  matrix.     Using 
Cramer's  rule  to  solve  the  above  linear  system  and  exploiting  the  prop- 
erty that  the  determinant  of  a  Vandermode  matrLx  generated  from  U\, . . . ,  u\f 
(denoted  by  V'(ui, .  . . ,  u.\/))  is  given  by  ni<j("i  ~  ";)•  ^'•^  obtain  that 

_    Wjllis) Ir-\{s).^.Ir  +  \{s) XSfjs))    _    -pj  Tfefs) 

''    ""  \-l  r.i  c\  T-w^cU  ~    li 


Therefore, 


Having  found  explicit  solutions  for  the  remaining  unknowns  Dr  we  have  explicit 
expressions  for  the  t„(5).  n  >  0-  The  proof  of  Theorem  1  is  now  complete. □ 
As  an  additianal  check  of  the  algebra  we  compute  the  generating  function 

Summing  up  (7)  and  (8)  we  obtain  that 

^M      (-l)-^3;^(0)      1-;         3,Jxr(s'\)-\      r-,.\f  ri,(5) 

For  r  =  1  we  obtain  that  ^{s.  1)  =  -,  which  is  the  condition  that  the  probabilities 
sum  to  one  in  the  transform  domain.  Another  interesting  point  is  the  fact  that 
the  generating  function  'l'(s.r)  is  symmetric  with  respect  to  the  roots  Xr{s).  This 
observation  is  nontrivial  and  it  is  established  by  using  the  Lagrange  interpolation 
formula  and  the  Chinese  remainder  theorem.  Theorem  I  was  proved  under  the 
assumption  that  the  roots  Xr{s)  are  distinct  and  therefore  all  the  formulae  are  valid 
only  in  this  case.  If.  however,  there  are  multiple  roots  (say  j:,(s)  =:  Xj{s)).  the 
resulting  formulae  are  simply  the  limit  of  the  formulae  given  here  as  x,(s)  — '  X]{s). 
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4      Busy  period  analysis 

Ramaswami  [12]  has  characterized  the  busy  period  of  an  G/Ph/\  queue  using  the 
matrix  geometric  approach.  In  this  section  we  simplify  his  result  considerably  and 
succeed  in  deriving  a  very  simple  formula  for  the  Laplace  transform  of  the  busy 
period  distribution  which  is  very  suitable  for  numerical  inversion  in  the  time  domain 
as  we  show  in  the  next  section. 

Let  Bp  be  the  random  variable  of  the  busy  period.  We  also  define  a  new  random 
variable  A,_,,  which  we  call  the  queue  build  up  time.  This  random  variable  plays 
a  critical  role  in  the  analysis.  A,,j  is  the  time  between  two  arrival  epochs  with 
the  following  properties.  Immediately  after  the  initial  arrival  epoch  there  are  n 
customers  in  the  system  and  the  STC  is  in  state  t,  while  immediately  after  the  final 
arrival  epoch  (there  may  be  other  arrivals  in  between)  there  are  n  +  1  customers  in 
the  system  and  the  STC  is  in  state  j.  In  addition,  throughout  the  time  A,,j  the 
number  of  customers  never  decreased  below  n.  Notice  that  A,,j  is  independent  of 
n. 

Let  R{t)  be  a  matrix  whose  i,j  element  is  the  pdf  R,.j{t)  =  ^Pr[Aij  <  t]  of 
the  queue  build  up  time.  Let  r,,j(s)  =  E[e~^^'-i]  be  the  transform  of  Ri,j{t)  and 
r(s)  be  the  transform  of  R{t).  Finally  let  ^'(s)  be  a  row  eigenvector  of  r(s)  with  an 
eigenvalue  u(s).  We  can  now  state  and  prove  our  basic  result. 


Theorem  2    The  Laplace  transform  cr[s)  of  the  busy  period  Bp  ts  given  by 

a(s)  =  E[e-'^n  =  1  -  (1  -  /?i(^))J;^;^'^^^;^,,,  (12) 

where  ir{s)  are  the  M  roots  of  the  polynomial  equation  (6). 


Remark:  (12)  is  a  simple  symmetric  function  of  Xr(s)  and  therefore  the  distinctness 

assumption  is  no  longer  necessary. 

Proof 
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By  focusing  on  the  last  customer  arrived  in  the  busy  period  we  write  the  busy  period 
dynamics. 

L  E  M".'(<)  *  {Mt)  *  6l""'V)  *  (  E  b\(t)Qr^l^]  r  a,{t)dt\  ,  (13) 

11=1  k=l  I  \r  =  l  /      •^'  J 

where  the  symbol  '*'  indicates  convolution,  the  first  term  corresponds  to  the  case 
in  which  the  last  customer  in  the  busy  period  was  the  only  one  in  the  busy  period 
and  in  the  second  term  (the  double  summation)  we  condition  on  the  number  of 
customers  n  and  the  state  k  of  the  STC  this  customer  has  found  when  he  entered. 
In  this  case,  the  busy  period  is  the  sum  of  two  independent  random  variables:  The 
time  to  build  n  customers  in  the  queue  (a  convolution  of  n  build  up  times  Ajjt) 
and  the  time  to  empty  the  system,  since  he  is  the  last  arriving  customer  in  the  busy 
period. 

Our  strategy  is  first  to  find  R{t)  and  then  using  (13)  to  find  the  transform  of 
the  busy  period.  We  are  thus  naturally  led  to  the  dynamics  of  the  queue  build  up 
times.  Similarly  to  the  first  passage  time  analysis  as  in  (Keilson  and  Zachmann  [8]), 
we  write  down  dynamics  of  the  queue  build  up  time  in  matrLx  form  by  considering 
the  last  arrival  during  A,,j 


mt)  = 


b\{t) 


b-^'{t) 


b\'f{t) 


a,(n 


r    /    ...,^    \ 


+  £/?("'(<)*< 


n=l 


bi{t) 

*  b[''-'\t)  *  (^  b\{t)     •••     b-l'{t)  )]ai(0  >  = 

=  S(<)ai(t)  +  ^  «(")(t)  *  [F„(Oai(0], 

n=l 

where  S(t),F„(t)  are  the  upper  diagonal  matrix  appearing  as  the  first  term  in  the 
sum  and  Fn{t)  is  the  matrix  composed  of  the  three  convolutions.    By  taking  the 
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Laplace  transform  of  the  above  matrix  equation  and  multiply  both  sides  from  the 
left  with  the  eigenvector  of  r(s)  ^'(s),  we  obtain: 

u{s)Cis)  =  ^{s)C{[B{t)  +  f;  u"(s)Fn(t)]ai(t)}- 


(14) 


n  =  l 


But, 


{Is  +  Bo  +  uBO~'  =  (/s  +  fio)-'  + 


1-U^l(5) 


3,{s)3,  is) 


Pm{s)I3i  is) 


since  for  every  pair  of  matrices  C  of  full   rank  and   D  of  rank   I,   (C  +  D)    ^    = 
^~^  ~  i  +  ir/c-'D'i  •   ^y  expressing  this  in  real  time  we  obtain 

b\{t)     ■■■     b'l'it) 


,-(Bo  +  uB,)( 


^   bi{t)   ^ 


n  =  l 


*b[''-'\t)*(^b\{t)     ...     b-l'{t)) 


\  b.Mit)  J 


0       ■••    b',',(t) 
As  a  result,  (14)  becomes 

u{s)^{s)  =  ^''(s)£{e-'S°  +  "(''5>)'a,(0}  =  ^{s)a,{Is  +  Bo  +  u(s)5i).  (15) 

Therefore,  since  ai(s)  is  a  rational  function,  ^(s)  must  be  a  row  eigenvector  of 
{Bo  +  u{s)Bi).  If  —:{s)  is  the  corresponding  eigenvalue,  following  the  same  technique 
as  in  claim  1  we  have  that 

ti(s)/?,(r(s))=l.  (16) 

Since  Qri(s)  is  a  rational  function  of  s,  we  get  from  (15)  that 

u(s)  =  a,(s--'(s)).  (17) 

Comparing  (16)  and  (17)  we  observe  that  :(s)  satisfies  exactly  the  same  equations 
as  x(s)  (equations  (6))  and  therefore  :{s)  =  x(s),  u(s)  =  w(s)  and  ^(s)  —  3\{x{s)). 
Having  characterized  the  eigenvalues  and  eigenvectors  of  T[s)  we  can  spectrum 
decompose  it  under  the  distinctness  assumption: 


M 


(X{s)r  =  X;(a,(s  -  x.(s)))"0r(s)/?i  (x.(s)), 


(18) 


r=l 
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where 


3{s)  = 


3i  (r,(5)) 


<^i(s)     •••     0\t{s) 


01  i^.sfis)) 

After  the  characterization  of  the  r(s)  we  take  the  Laplace  transform  of  (13). 
Using  (18)  and  after  similar  manipulations  with  the  analysis  of  the  queue  build  up 
time,  we  obtain 


r  =  l 


=       ",(^(5)) 


-1 


/     ,3i(ri(5))-l     \ 


/?l(xM(a))-l 


Since  for  any  non  singular  matrix  .4  we  know  that  f  .4    ^y    =    1    —      ^j  ,//f    .  we 
find  that 


<t(s)=  1- 


1 


det(^(5)) 


det(/3(s)- 


,-f  3i(xi(s))-l 
^1       »-r,l5) 


^  3i(rAf(0)-l 

1        S-X.Vf(s) 


)• 


But  for  all  r=  1,...,A/ 

^3x{zr{s))-\ 

1 

^'        S-Ti^s) 

S  -  Xr{s) 

1 

S  -  Xr^s) 

1 

S  -  Xr{s) 

1 

S  -  Xr(s) 


{e-'i5,(x.(5))-e-l} 


e-;{(/r.(5)+5o)-n-5i)-/} 
e-'i(/x.(s)  +  flo)~'{-5,  -  (/x,(s)  +  So)} 


Ji(r,(s)){/(s  -  Xr{s))  -{Is  +  Bi  +  Bo)} 


Hence, 


<T{S)       =       1- 


=     Ji(x.(s)) 3,{xris)){Is  +  5i  +  5o). 

S  -  Xr[s) 


det( J(5)  -  { J(5)  -  diagonal(;3J^,  •  ■  • ,  jZ7;;iJ^)g{s)iIs  +  B,  +  Bq)}! 

det(J(s)) 
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=     1  -det(diagonal( r-T'--- T":))  ^«'t(/5  +  Bi  +  Bq) 

s-xi(s)  s-x\f[s) 

Denominator  Ji(s)  —  Numerator  Ji{s) 


^2i';'(l-Q^)/^0^•-•^l^^)) 


Although  the  analysis  used  some  rather  heavy  machinery  from  linear  algebra, 
it  used  direct  probabilistic  arguments  by  considering  the  dynamics  of  the  system. 
The  reward  of  this  analysis  is  a  very  simple  expression  for  the  transform  of  the  busy 
period  distribution,  which  as  we  show  in  Section  6  offers  very  important  computa- 
tional advantages.  In  addition,  it  is  not  hard  to  compute  in  closed  form  moments  of 
the  busy  period  distribution  by  repeated  differentiation  of  the  Laplace  transform. 

5      The  waiting  time  distribution  under  FCFS 

In  this  section  we  derive  an  expression  for  the  conditional  waiting  time  pdf  given  the 
arriving  time.  Our  strategy  for  the  analysis  is  the  following;  we  first  find  the  distri- 
bution of  the  number  of  customers  in  the  system  at  an  arrival  epoch;  conditioned 
on  the  number  of  customers  found  upon  arrival,  we  then  find  the  waiting  time  pdf 
and  finally  we  find  the  (unconditioned)  waiting  time  pdf. 

In  order  to  obtain  a  closed  form  expression  in  the  transform  domain,  we  make  the 
assumption  that  the  system  is  initially  empty  and  futhermore,  the  initial  probability 
distribution  has  a  very  special  form. 

Assumption  1    We  assume  that  the  initial  probability  vector  Po{Q)  =  Aai(O). 

In  principle,  this  assum[)tion  is  not  necessary  if  we  take  a  pure  numerical  approach 
for  the  solution.  Without  this  assumption,  however,  it  is  not  possible  to  obtain  a 
closed  form  expression  both  in  real  time  and  in  the  transform  region.    In  the  next 
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theorem  we  prove  a  critical  consequence  of  assumption   1;  the  arrival  process  has 
already  reached  steady  state  from  the  beginning,  i.e., 

Po,it)  =  Pr[Rait)  =  (]  =  Pr[i?a(0)  =  i]  =  Po,.(0). 

Proposition  1  If  the  initial  condition  satisfies  assumption  J,  the  arrival  time  of  the 
first  customer  IS  the  forward  recurrance  time  of  the  tnterarrival  time,  i.e.  the  residual 
life  time  of  the  renewal  interval.  Furthermore,  Po(0)  is  the  stationary  solution  of 
the  Kolmogorov  equation: 

±P\t)  =  -P'{t){Ao  +  A,) 


P'[t)    1  =  1, 


that  describes  the  arrival  process. 


Proof 

The  transform  pdf  of  the  interarrival  time  T^  of  the  first  customer  is  given  by. 


a'{s)  =  E[e''^]  =  P^{0) 


I  ai(5)  ^ 


V  ^U«)  / 


^   ai(5)  ^ 


Since  ai(0)   =  e*i/lo      ^"^ 


=   {Is  +  .4o)~'(-.4ie-i),   we  find  after 


simple  algebraic  manipulations  the  transform  of  the  forward  recurrance  time  under 
assumption  1  as  follows: 

(x'{s)     =     Aa;(0)-(/s  +  .4o)-'(-.4ie-'i) 
=     \^^{A^)-\U^A^)-\-Axe^) 
=     Ae",!  ((^o)-'  -  {Is  +  .40)-')  (-.4ie-i) 

=     -(l-ai(s)). 
s 
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To  obtain   the  stationary  distribution,   we  observe  that  (.4o  +  --^i)!    =   0  and 
because  of  the  structure  of  .4i ,  we  know  that  .4i  1  =  i4ie*i.  Thus  we  find 

l^-Ao^Aiex,  (19) 

or  equivalently 

We  are  now  ready  to  prove  that  the  stationary  probability  vector  is  proportional  to 

ai(0);  for  this  it  suffices  to  show  q'j(0)(.4o  +  .4i)  =  0'.    Since  q'i(O)  =  ^^Aq^,  we 

have 

q'i(0)(.4o+     .4i) 

=     e-*i.4o'(.4o  +  ^i) 

=     e*i  +  e*i  /Iq  '  ^  1 

=    0' 

We  complete  the  proof  by  showing  j  =  <^i(0)  ■  1-    Utilizing  (19)  and  the  definition 
a-[{s)  —  — e*i(/s  +  .4o)~'.4ie'i  we  obtain 

a;(o)-   r 

=     -a',{0)A-' Aiei 

=     -\[ms-,oe'iAQ\ls  +  Ao)~'^Aie'i 
=     lim,_oie-\((/s  +  .4o)-i-.4o-').4,e-, 

_     i 

Therefore,  we  have  proved  that  the  ergodic  solution  to  the  above  Kolmogorov  equa- 
tion is  Aq,(0)  =  Po(0).  □ 

A  corollary  of  the  theorem  is  that  the  expression  for  Dr  in  theorem  1  further 


22 


simplifies  to 


°-  =  7  ii"(x,(.))  "Wn, .(„_,.(„■  (20) 

since 

^  A 

I]  Poa(OW(s)  =  a'(s)  =  -  (1  -  a,(s)) . 

We  will  next  find  the  distribution  of  the  number  of  customers  in  the  system 
seen  by  an  arriving  customer.  We  define  the  event  .4.40  —  Arrival  about  to  occur  in 
(r,  r  +  dr)  and  the  pre-arrival  probabilities:  P^ti''')  =  Pr['^'(^)  =  "■  ^a('')  =  i|/1^0]- 
In  the  following  proposition  we  find  the  pre-arrival  probabilities. 

Proposition  2    Under  assumption  1  the  vector  of  the  pre-arrival  -probabilities  is 


and  its  Laplace  transform  is 

^n{s)     =     {^D.Ji(x.(s))(q,(5-x.(5)))"  {n>  1} 


1  '^ 


^.=1 


^0{S)       =       TY,Dr3]{Xr{s)). 


^.1 


Proof 


_  Pr[uf^i(:V(r)  =  nDRAr)  =  inR^jr)  =  /)  fl  .4.40]  _ 
Pr[ujL,i?a(r)  =  /n.4.40] 

ELi  Pr[.4.4C'|.V(r)  ^  n  H  /?,(r)  =  i  H  fl,(r)  =  /]  Pr[.V(r)  =  n  H  /il.fr)  =  i  n  /i:,(r)  =  /] 

ELi  Pr[.4.40|/?a(r)  =  /]  Pr[fla(r)  =  /] 
But 

Pr[.4.40|.V(r)  =  nnRs{r)  =  in  Ra[r)  =  /]  =  Pr[AAO\Ra{T)  =  I]  =  Xipidr 

and 

Pr[/t,(r)  =  /]  =  Po,/(0)  =  Acti(O), 
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from  proposition  1.  It  is  at  this  point  that  assumption  1  becomes  critical.  Without 
this  assumption  Pr[/?a(^)  =  ^  would  be  a  function  of  t.  wiuie  under  assumption  1 
it  is  independent  of  r.  and  therefore  it  would  not  be  possible  to  find  a  closed  form 
formula  for  the  transform  of  the  pre-arrival  distribution.  Therefore,  we  obtain 

since  ^i-i  ^ipio:[{0)  =  fti(O)  =  1.  Therefore,  using  vector-tensor  notation  we  have; 

pr(^)  =  {^;(o{(-.4,e-,)®/}. 

In  the  transform  region,  using 

r  =  l 

we  obtain  (the  derivation  for  ttq^s)  is  similar) 


M 


^n(-0     =      -^Dr/?i(r.(5))(«i(5-r.(.s))r   {n>  1}, 


r=l 


^o"(«)    =    ^l]a/?J(r.(s)).n 


r=l 


We  are  now  ready  to  prove  the  central  theorem  of  this  section. 

Theorem  3    Under  assumption  1  (he  Laplace  transform,  of  the  waiting  time  distri- 
bution under  FCFS  is 


f 

Jo 


1 


M 


M 


e-  Prmr)  <  t]dr  =  i  +  ^  ^"X'^u'  I  11 


■rj5) 


\ 


S  ~;     S(3-^'{lr{s))      \t=\    ^r{s)-  Xk{s) 


.^r{s)t 


Proof 


Given  there  are  exactly   n  customers  in   the  system  including  the  customer  just 
arrived  and  the  STC  is  in  stage  i,  then  waiting  time  c.d.f.  is: 

/J  6.(0  *  6l"-''(0  ♦  Ejl,  b\{t)q,fi,dt     when  n  >  3 

/o  T.'jL\  bi{t)qjfijdt  when  n  =  2 

,  U{t)  when  n  =  1 
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where  U{t)  is  unit  step  function.  By  conditioning  on  the  state  the  arriving  customer 
found  the  system,  and  using  the  expressions  for  the  pre-arrival  probabilities  from 
Proposition  2  we  obtain: 


f^     e-'^PT[lV{T)<t]dT 


I 


+   Er=i(ai(s-x.(s)) 


in+l 


ti(0 


h\(t) 


*'''r''(0*(    6}(0 


^   QlMi   ^ 


V  ^.\K<)  / 


di\ 


21) 


We  have  observed  in  the  analysis  of  busy  period  from  the  previous  section  that 
h\{t)     ...     h\\i) 


-{So  +  uS])(  _ 


0     ■■■   fe:l{(0 


^   6i(0   ^ 


+E"" 


n  =  l 


*fc';~''(0*(   61(0 


,.\/ 


h\'[t) 


\  bxtit)  I 


Substituting  this  to  (21)  we  obtain 

/'     e-''PT[W{r)  <t]dr  = 
Jo 

-  2  Dr3[{xr[s))  i^ed'it)  +  j^  a,{s  -  x.(s))e-(^°+"'(^-^^(^''^'"(-5ie-i)(i< 

Since  3\{^Xr{s))  is  a  row  eigenvector  of  [Bq-\- ci\{s  —  Zr[s))B\ )  with  eigenvalue  — Xr(s) 
and  ai(s  —  Xr(s)),'^i(-Cr(s))(  — Sie*!)  =  1,  we  obtain 

/o^     e-^^  Pr[ir(r)  <  <](ir 

=     \  YM.  DrUit)  (^l(r.(5))e-  +  ^i^) 

=     }  E.^I.  aC'-(t)  {3[{xrisW  -  ^^^^f^).-:  +  ^) 

=  iE.'ii  ar(o  (-^,J;(x.(s))floe'i  +  ^) . 
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Substituting  expression  (20)  for  D^  and  using  (11)  we  finally  obtain 


Vf 


1    ,     -    (-l)V'J,vr(0)       -  r,{s) 


e'^'<'''.n 


A  corollary  from  the  previous  theorenn  is  that 

$(s,-.')     =      /"  /"  e-'^-^'-^  Pr[W{T)  <  t]dtdT 
Jo     Jo  en 

Note  that  this  expression  is  a  synnmetric  function  of  the  j;r(s)'s.  In  addition  the 
quantity  linij^o  s<t>(s,u;)  is  the  solution  of  the  steady  state  Lindley  equation  for  the 
G//G/1  queue,  i.e.  the  transform  of  the  steady  state  waiting  time  distribution. 

6      Numerical  Results 

In  the  previous  sections  we  have  derived  explicit  expressions  for  the  Laplace  trans- 
forms of  the  queue  length,  the  waiting  time  and  the  busy  period  distributions.  In  this 
section  we  will  remove  the  "Laplacian  curtain",  by  numerically  inverting  the  Laplace 
transforms.  The  numerical  inversion  of  the  Laplace  transform  is  a  well  studied  but 
not  completely  solved  problem  in  numerical  analysis.  In  fact,  Platzman,  Ammons, 
Bartholdi  [11]  show  that  the  problem  of  numerically  inverting  the  Laplace  transform 
of  a  probability  distribution  is  #P  complete,  that  is  a  hard  computational  problem. 
Our  overall  algorithm  is  written  using  the  software  package  of  Mathematica, 
developed  by  Wolfram  [13]  and  works  2is  follows.  We  first  compute  the  roots  of  the 
polynomial  equations  (6)  for  selected  s  values.  For  this  purpose  we  use  the  build  in 
functions  of  Mathematica  to  find  all  the  roots  of  (6).  We  then  use  the  algorithms 
of  [11]  and  [7]  to  compute  the  inverse  Laplace  transform  of  the  distributions  under 
study.  We  used  two  algorithms  to  invert  numerically  the  Laplace  transforms: 

1.   The  algorithm  of  Platzman  et.  al.  [11] 

This  algorithm  works  for  distributions  that  are  defined  over  finite  regions.  We 
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used  this  algorithm  combined  with  fast  Fourier  transform  for  the  inversion  of 
the  busy  period  distribution.  Although  the  busy  period  takes  values  in  the 
region  (0,oc),  in  order  to  apply  the  algorithm  we  used  the  region  (0,  E[Bp]  + 
3Var[Bp])  as  the  region  on  which  the  busy  period  is  different  from  0.  For 
details  of  this  algorithm  the  reader  is  referred  to  [11], 

2.  The  algorithm  by  Hosono  [7] 

Hosono  [7]  proposed  an  algorithm  for  inverting  Laplace  transforms  which  is 
quite  robust  and  accurate.  We  used  this  algorithm  for  numerically  inverting 
the  busy  period,  the  waiting  time  and  the  queue  length  distributions.  The 
algorithm  is  not  well-known  in  the  western  literature  since  it  is  primarily  pub- 
lished in  Japanese.  We  found  however  that  it  is  a  very  robust  algorithm.  It 
satisfies  all  the  necessary  conditions  which  an  ideal  fast  algorithm  should  sat- 
isfy, namely:  It  is  easy  to  program  and  control  the  error.  Moreover,  it  has 
small  memory  requirements,  short  computational  times  and  can  be  used  for  a 
wide  variety  of  problems.  We  briefly  introduce  the  algorithm  below. 

Let  <i>{s)  be  the  input  function.  Let  f{t)  be  the  inverse  Laplace  transform  of 
(t>{s).  We  choose  a  precision  p  {significant  digits}  so  that  error  of  numerical 
inversion  is  less  than  10~^'''^|/(<)|.  Let 

The  algorithm  works  as  follows: 


(a)   Find  it  so  that 


E  :^- 


r=0 


< 


.2 


(b)   Evaluate  f{t): 


k-\  p-i 

n=l  r=0 
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Hosono  claims  that  this  algorithm  works  when  f{t)  is  sufficiently  smooth'. 

Comments  on  the  numerical  results 

All  the  computation  was  done  in  a  Macintosh  II,  and  all  the  program  is  written 
in  Mathematica.  For  computing  the  transient  queue  length  and  the  waiting  time 
distribution,  we  assumed  that  the  arrival  time  of  the  first  customer  is  the  forward 
recurrence  time  of  the  interarrival  distribution. 

From  our  preliminary  experience  we  can  say  that  the  algorithms,  in  particular  the 
algorithm  by  Hosono,  are  robust  and  run  very  fast.  The  largest  example  we  ran  was 
an  MGE2o/MGE2o/\.  which  was  solved  by  the  algorithms  without  any  difficulty. 
Unfortunately,  we  did  not  have  any  other  numerical  results  to  compare  with  except 
the  ones  for  the  A//A//I  queue  queue,  whose  solution  is  known  explicititly  in  terms 
of  modified  Bessel  functions  (see  Gross  and  Harris  [6],  p. 143).  In  section  61  we 
compare  the  results  of  this  algorithm  to  the  exact  results.  As  an  illustration  of  the 
stability  and  robustness  of  the  algorithms  we  present  in  section  6.2  an  example  of 
the  MCE3/MCE2/I  queue. 

6.1      The  A//A//1  queue 

We  start  our  examples  with  an  A//A//1  queue  with  traffic  intensity  p  =  0.75.  The 
interarrival  rate  is  A  =  1  and  service  completion  rate  /i  =  4/3.  In  order  to  compare 
the  accuracy  of  the  two  algorithms  with  the  known  solution  we  computed  in  table  I 
the  CDF  of  the  busy  period  (The  mean  is  E[Bp]  =  3  and  the  coefficient  of  variation 


Fn  needs  to  satisfy  the  following  conditions: 

•  for  sufficiently  large  n,  1/2  <  \F„+JFr,\  <  I. 

•  when  n  —  oo,  fn,  Af„,  A^fn,  •  •  ■  converge  monotonically  to  0, 

where  A"^  denotes  the  r  th  difference. 

It  can  be  shown  that  the  violation  of  these  conditions  results  in  the  Gibbs  phenomenon,  which  only 

appears  at  points  of  discontinuity  of  f{t). 


28 


is  Cg   =   7.)  The  algorithm  by   Hosono  gives  identical   results  with  the  the  exact 
solution. 

In  figure  2  we  plot  the  first  and  second  moments  of  the  waiting  time  as  a  funtion 
of  ^  In  figure  3  we  plot  the  waiting  time  distribution  as  a  funtion  oft.  In  figures  4, 
5  we  plot  the  first  and  second  moments  and  the  distribution  of  the  queue  length  as 
a  funtion  of  t. 

6.2     The  MGE3/MGE2/].  queue 

Merely  as  an  illustration  of  the  algorithms  we  chose  a  MCE3/MGE2/I  queue  with 

the  following  distributions:  The  interarrival  distribution  has  Laplace  transform: 

2                             2         4  2         4         6 

Qi(s)  =  0.5- +  0.5  X  0.3- +  0.5  x  0.3; 


2+s  ■    2+s4+s  2+s4+s6+s 

with  mean  E[T]  =  0.683333  and  coefficient  of  variation  is  Cj-  —  0.70. 
The  service  distribution  has  Laplace  transform: 

5  5        3 

/3ifs)  =  0.8- +  0.2- 


5+s  5+s3+s 

with  mean  ^^[A']  =  0.266667  and  coefficient  of  variation  is  C^   =  1.125;  thus  the 
traffic  intensity  is  p  =  039. 

By  differentiating  the  transform  of  the  busy  period  distribution  we  found  that 
the  mean  of  the  busy  period  is  E[Bp]  =  0. 398444  and  the  coefficient  of  variation 
is  Cg  =  2.327.  In  figure  6  we  plot  the  busy  period  CDF,  in  figure  7  we  plot  the 
first  and  second  moments  of  the  waiting  time  as  a  funtion  of  t.  In  figure  8  we  plot 
the  waiting  time  distribution  as  a  funtion  of  ^  In  figures  9,  10  we  plot  the  first  and 
second  moments  and  the  distribution  of  the  queue  length  as  a  funtion  oft. 

7      Concluding  Remarks 

In  this  paper  we  attempted  to  demonstrate  the  power  of  root  finding  techniques 
for  problems  which  were  considered  intractable  like  the  queue  length  and  waiting 
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time  distribution  in  the  transient  domain  and  the  busy  period  distribution  for  the 
MCEil MGE\i/\  queue.  Using  direct  probabilistic  arguments  combined  with  tech- 
niques from  linear  and  tensor  algebra,  we  succeeded  in  deriving  explicit  expressions 
for  the  Laplace  transforms  of  the  distributions  under  study. 

Algorithmically  our  approach  offers  a  method  for  finding  these  distributions  in  the 
time  domain  through  the  numerical  inversion  of  the  Laplace  transforms.  Our  ex- 
periments with  the  method  are  very  encouraging  since  our  experience  with  the  al- 
gorithms we  used  was  very  positive,  since  they  are  very  fast,  robust  and  very  easily 
programmable. 
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t 

Exact 

PA.B 

Hosono 

0. 

0. 

0. 

0. 

0.94104 

0.565744 

0.529196 

0.565744 

2.07421 

0.734767 

0.732878 

0.734767 

3.20737 

0.806977 

0.811223 

0.806977 

4.34054 

0.848456 

0.854939 

0.848456 

5.47371 

0-875836 

0.883652 

0.875836 

6.60688 

0.875836 

0.904243 

0.875836 

7.74004 

0.9102 

0.919855 

0.9102 

8.87321 

0.921766 

0.932158 

0.921766 

10.0064 

0.931074 

0.942134 

0.931074 

11.1395 

0.938727 

0.9504 

0.938727 

12.2727 

0.945129 

0.957369 

0.945129 

13.4059 

0.950558 

0.963326 

0.950558 

14.539 

0.955216 

0.968478 

0.955216 

15.6722 

0.959252 

0.972977 

0.959252 

16.8054 

0.962779 

0.976937 

0.962779 

17.9386 

0.965883 

0.980449 

0.965883 

19.0717 

0.968631 

0.983583 

0.968631 

20.2049 

0.97108 

0.986393 

0.97108 

21.3381 

0.973271 

0.988927 

0.973271 

22.4712 

0.975241 

0.991219 

0.975241 

23.6044 

0.977019 

0.993303 

0.977019 

24.7376 

0.97863 

0995202 

0.97863 

25.8707 

0.980094 

0.996939 

0.980094 

29.2702 

0.983766 

1. 

0.983766 

Table  I:  The  A//A//1  busy  period  CDF. 
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Figure  2:  The  first  and  second  moments  of  the  waiting  time  of  an  A//M/1  queue  as 
a  function  of  time 
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Figure  3;  The  waiting  time  distribution  of  an  M/M /\  queue  as  a  function  of  time 
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Figure  4:   The  first  and  second  moments  of  the  queue  length  of  an  M/M/l  queue 
as  a  function  of  time 
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Figure  5:  The  queue  length  distribution  of  an  M/M/\  queue  as  a  function  of  time 
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Figure  6:  The  busy  period  CDF  of  an  MCE3/MCE2/I  queue 
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Figure  7:  The  first  and  second  monnents  of  the  waiting  time  of  an  MGE3/MGE2/I 
queue  as  a  function  of  time 
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Figure  8:  The  waiting  time  distribution  of  an  MCE3/MCE2/I  queue  as  a  function 
of  time 
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Figure  9:  The  first  and  second  moments  of  the  queue  length  of  an  A/G£'3/A/G£'2/l 
queue  as  a  function  of  time 
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Figure  10:  The  queue  length  distribution  of  an  MCE2/MCE2/I  queue  as  a  function 
of  time 
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